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The isospin-2 tttt system provides a useful testing ground for determining elastic hadron scat- 
tering parameters from finite- volume spectra obtained using lattice QCD computations. A reliable 
determination of the excited state spectrum of two pions in a cubic box follows from variational 
analysis of correlator matrices constructed using a large basis of operators. A general operator 
construction is presented which respects the symmetries of a multi-hadron system in flight. This is 
applied to the case of tttt and allows for the determination of the scattering phase-shifts at a large 
number of kinematic points, in both S-wave and D-wave, within the elastic region. The technique 
is demonstrated with a calculation at a pion mass of 396 MeV, where the elastic scattering is found 
to be well described by a scattering length parameterisation. 

PACS numbers: 14.40.Be, 12.38.Gc 



I. INTRODUCTION 

QCD is the accepted underlying theory of the strong 
interactions, and the properties of the spectrum and in- 
teractions of hadrons should be calculable from it using 
a suitable regularization scheme for the quark and gluon 
fields. A particularly convenient approach is to consider 
the theory on a finite lattice of space-time points so as to 
admit a numerical method of solution. While significant 
progress has been made recently in determining the single 
particle spectrum of hadrons, describing the resonances 
seen in scattering experiments in terms of eigenstates of 
QCD has remained a challenge to lattice calculations. 
Direct access to the matrix elements related to decays is 
missing in the Euclidean formulations of lattice QCD. In 
principle, the relevant hadronic matrix elements can be 
inferred indirectly through a detailed study of the spec- 
trum in a finite- volume lattice box [2 [5]. Within this 
approach, one can map the discrete spectrum of eigen- 
states of the finite volume theory to the infinite volume 
scattering parameters, and if present, observe resonant 
behavior. 

Crucial to this approach is the high-precision determi- 
nation of multiple excited eigenstate energies with a given 
quantum number. Determination of the discrete spec- 
trum of finite- volume eigenstates follows from analysis of 
the time-dependence of two-point correlation functions 
featuring operators of the desired quantum numbers con- 
structed from quark and gluon fields. For creation and 
annihilation at time and t respectively we have 

a,(i) = (o|o,(i)o](o)|o). 

Inserting a complete set of eigenstates of the Hamilto- 



nian, this correlator has a spectral decomposition 
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where the sum is over all states that have the same quan- 
tum numbers as the interpolating operators Oi , Oj . Note 
that in a finite volume, this yields a discrete set of en- 
ergies. En- It is these finite- volume energies that are 
related to infinite volume scattering amplitudes through 
the Liischer method [5]. 

A relatively straightforward sector in which to study 
hadron scattering in finite-volume is tttt in isospin-2. At 
low energies, this channel is experimentally observed to 
be non-resonant in low partial- waves [3 6J and this lack 
of resonances ensures a slow variation of phase-shifts 
with energy. This makes the problem of determining the 
phase-shift as a function of energy somewhat easier. A 
difficulty of this choice of channel is that the interaction 
between pions in isospin 2 is weak so that the discrete 
energies in finite-volume are shifted relatively little from 
the values relevant for non-interacting pions. This will 
require us to make precision measurements of the en- 
ergy spectrum in order to resolve the differences. Within 
the field-theory, the correlators for this channel do not 
contain any annihilation contributions; the only Wick 
contractions featuring are those in which the four quark 
fields in the creation operator (at t = 0) propagate in 
time to the annihilation operator (at t). The absence of 
quark propagation from i to i reduces the computational 
overhead for the calculation. 

In a previous publication [7] we presented the first lat- 
tice QCD study of the energy-dependence of S and D- 
wave TTTT scattering in isospin-2. We limited ourselves to 
the TTTT system overall at rest and found only a handful of 
points below the An inelastic threshold on the lattice vol- 
umes considered. In this paper we will also consider the 
TTTT system "in-flight" , that is with an overall momentum 
(satisfying the periodic boundary conditions of the finite 



cubic lattice). This allows us to determine the phase- 
shifts at a larger number of discrete energies below the 
47r inelastic threshold and to map out the energy depen- 
dence of the scattering in more detail. The price to be 
paid is that the relevant symmetry group in the lattice 
calculation is significantly reduced. At rest the lattice 
has a cubic symmetry whose irreducible representations 
("irreps") contain multiple angular momenta, e.g. the 
"scalar" representation, Af contains as well as £ = 0, 
also i = A and higher. In-flight, with two pious having 
total momentum, P, the symmetry is restricted to rota- 
tions and reflections which leave the cubic lattice and the 
axis defined by P invariant. The irreps of this symmetry 
group are even less sparse in tttt scattering angular mo- 
mentum; the "scalar" representations typically contain 
^ = 0,2,4.... In this work we will consider the effect 
these higher partial waves have on the determination of 
scattering phase-shifts for the lowest £ values. 

In [7], we used only the simplest V'TsV' interpolators 
in construction of tttt correlators. Single-pion correla- 
tors constructed with these operators are saturated by 
the ground state only at rather large times, and similarly 
the TTTT correlators receive significant contributions from 
excited tt* states. The need to consider correlators at 
large times increases the degree to which we feel the sys- 
tematic effect of the finite temporal extent of the lattice 
(T) . Limited account was taken of these effects in [7] . In 
this paper we take steps to address finite-T effects, firstly 
by using "optimised" pion operators which are saturated 
by the ground state pion at earlier times, and secondly 
by explicitly attempting to remove the leading effects of 
finite-T from the measured correlators. While these ef- 
fects are small in absolute terms, determination of the 
rather weak I ~ 2 interaction relies upon precise mea- 
surement of small energy shifts, and as such it is impor- 
tant to account for even small systematic effects. 

Our approach to determining the finite-volume spec- 
trum is to use a large basis of operators in each symme- 
try channel with which we form a matrix of correlation 
functions having all relevant operators at the source and 
sink. This matrix can be analysed variationallvjHHTU] . 
extracting a spectrum of energy eigenstates which are 
orthogonal in the space of operators used. This orthogo- 
nality is particularly useful in cases where levels are close 
to degenerate and to extract relatively high-lying states 
whose contribution to any single correlation function may 
be small relative to the ground state. The excited single- 
hadron spectrum of isovector and isoscalar mesons [II]- 
[15] and baryons [TU [IS] has been extracted with some 
success using this procedure. In the present case we re- 
quire a basis of operators capable of interpolating a pair 
of pious from the vacuum, constructed to transform ir- 
reducibly in the relevant symmetry group. The fact that 
J = 2 is expected to have only relatively weak inter-pion 
interaction strength suggests a natural basis might be 
one resembling pairs of non-interacting pions, i.e. pions 
of definite momentum. 



In general our tttt creation operators have the form 



/ \[fel,fc2]t 

(TTTT) s, . 



E 

fcl,fc2 
fcl+fe2=-P 



C(P, K IX- kl-MT^\kx) T:\k2) 



here C are the Clebsch-Gordon coefficients for combining 
the two pion operators of definite momentum fci, hi so 
that the operator overall transforms in the irrep A of 
the relevant symmetry group for total momentum P — 
fci + fc2. This involves summing over multiple values of 
momenta fci, ki with the same magnitudes, \k\\, \k2\ and 
related by allowed lattice rotations. The basis is built 
up out of different magnitudes of pion momenta that can 
sum to give the same P. Much greater detail will be 
presented later in this paper. 

Using this basis we compute correlators within various 
irreps A for various P, 
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and for a fixed P, A, fi we perform variational analysis 
in a basis of operators labeled by [fci,fc2] leading to a 
finite-volume spectrum, En{P,A;L). This spectrum, de- 
termined in the rest frame of the lattice, corresponds to 
a discrete set of scattering momenta, Pcm j in the center- 
of-momentum frame. 

The finite-volume spectrum so obtained is related 
through the Liischer formalism [21 E] (as extended in [161- 
[T5] to the case of moving frames) to the phase-shifts, 
di{pcm), for elastic tttt scattering in partial waves of an- 
gular momentum, i. As discussed earlier, a given irrep 
A of momentum P, contains multiple angular momenta, 
£, and the formalism relates the finite-volume spectrum 
to the scattering amplitudes for all relevant £ though the 
following formula: 



det 



,2i5(pcm) 



U^^'^HPc.^] 



= 



(2) 
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Here U'^^''^)(pcm^) is a matrix of known functions and 
g2j(5(pcm) jg g^ diagonal matrix featuring the scattering 
phase-shifts {5i}. In both cases the rows and columns 
of the matrices are labelled by the angular momenta, £, 
relevant for the irrep (P, A). These matrices are for- 
mally infinite, but we may take advantage of the hierar- 
chy 6o ^ 62 ^ Si . . . relevant at low energies 
to reduce the effect of higher £ in Equation 

We will explore two methods to extract the phase 
shifts. The first method, similar to the one used in [7], ex- 
ploits the above hierarchy to determine the phase-shift in 
the lowest contributing partial wave and estimates a sys- 
tematic uncertainty from plausible variation of the higher 



^ near threshold, angular momentum conservation requires Si 

Pern 



partial wavesjj The second method parameterizes the 
momentum-dependence of the phase-shifts in £ = 0, 2 . . . 
using effective range expansions, then by performing a 
global fit which attempts to describe many finite- volume 
momentum points in many irreps, finds the values of the 
effective range expansion parameters. 

Computations were performed on anisotropic lattices 
with three dynamical flavors of Clover fermions 20j 2fj 
with spatial lattice spacing Ug ~ 0.f2fm, and a temporal 
lattice spacing approximately 3.5 times smaller, corre- 
sponding to a temporal scale a^^ ^ 5.6 GeV. This fine 
temporal lattice spacing has proven useful in determin- 
ing the spectrum of mesons and baryons, as well as the 
previous tttt I — 2 results [7]. In this work, results are 
presented for the light quark 1 = 2 spectrum at quark 
mass parameter atmi = —0.0840 and atms = —0.0743 
corresponding to a pion mass of 396 MeV, and at lattice 
sizes of 16^ X f28, 20^ x f 28 and 24^ x 128 with corre- 
sponding spatial extents L ^ 2fm, ^ 2.5 fm and ^ 3fm. 
Some details of the lattices and propagators used for cor- 
relation constructions are provided in Table |lj 

Recently, the NPLQCD collaboration [22] has deter- 
mined the £ = scattering phase-shift on the same en- 
sembles as used in this study, plus an additional larger 
lattice volume ~ 4fm. Their calculation is limited in 
scope by the fact that their approach does not project 
pion operators of definite relative momentum at the 
source. We will compare the results of the different 
approaches later in this paper. Other studies of tttt 
1 = 2 scattering in lattice QCD ([23-25J) have largely 
limited themselves to the threshold behavior of the scat- 
tering amplitude in S'-wave, as expressed by the scatter- 
ing length. 

Readers who are not concerned with the details of the 
calculation can skip to Section [IX| where the results for 
elastic scattering are presented. The remainder of the 
paper is structured as follows: 

Section |ll] outlines the construction of a basis of irre- 
ducible TTTT operators at rest and in-flight from products 



{L/usf X (T/at) 



of pion operators of definite momentum. Section III de- 



scribes the construction of correlators using the distil- 
lation framework. Section |IV| presents "optimised" sin- 
gle pion operators constructed as linear combinations of 
composite QCD operators with pion quantum numbers. 
Section |V] discusses the determination of the pion mass 
and anisotropy from measurements of the pion disper- 
sion relation. Section I VII considers the effects of the fi- 
nite temporal extent of the lattice on tttt correlators and 
presents mechanisms for reducing the role of these ef- 
fects in the determination of the discrete energy spec- 
trum. Section |VII| presents the finite- volume spectrum 
obtained on three volumes. Section IVIIII discusses the 
extraction of elastic scattering phase-shifts from finite- 
volume spectra using the Liischer formalism including a 



le-" X 128 
20^ X 128 
24^ X 128 



Nctss iVt,,„ A^v. 



479 12 



601 
553 



64 
128 



5 (P=0) 
3 (Pt^O) 

3 162 



TABLE I. The lattice ensembles and propagators used in this 
paper. The light and strange quark mass are atirii = —0.0840 
and atrUs = —0.0743 described in Ref. [21], corresponding 
to a pion mass of 396 MeV. The lattice size and number of 
configurations are listed, as well as the number of time-sources 
and the number of distillation vectors Vvccs (to be described 



in Sec III I featuring in the correlator construction. 



study of the possible effect of sub-leading partial waves 
within a toy model. Section |IX| presents our results for 
(5o,2 in the region of elastic scattering for 396 MeV pi- 
ous. Section |X] summarises our approach and results and 
suggests future applications of the methodology. 



II. OPERATOR CONSTRUCTION 

In order to calculate scattering amplitudes we must 
extract multi-hadron energy levels with high precision 
and so need interpolating operators that efficiently in- 
terpolate these multi-hadron states. To achieve this we 
consider operators constructed from the product of two 
hadron operators projected onto definite total momen- 
tum, P, and transforming as a definite irreducible rep- 
resentation of the appropriate symmetry group, lattice 
irrepnA, with irrep row, fi, 



Oa^(P) 



t 



y^ C{PAfi;kiAifii;k2A2H2) 



ki,k2 
k-i_+k2=P 



X OAi^i(fci) 



5'A2M2(fc2 



(3) 



Here Oai,^i (fci) and Oaj^^^ (^2) are hadron operators (for 
example, fermion bilinear operators), each projected onto 
definite momentum, irrep and irrep row. The Clebsch- 
Gordan coefficients, C, and the momenta appearing in 
the sum over ki and ^2 will be discussed later. 

A conventional infinite volume continuum analogue of 
this construction (for total momentum zero with p = 
fci = — ^2) would be 



O 



[SA 

J,M 



Ai A2 



J2 dp C{JiSM;pSi\i;~p 82X2) 



lmS2^2 



[o^'^'{p)y [o'^'^^-p)y , 



(4) 



we note that Ref. |19| has recently discussed a similar approach 



we use 'lattice irrep' to refer to the octahedral group irrep for 
a particle at rest and the irrep of the appropriate little group, 
discussed later, for a particle at non-zero momentum 



with 

C^{SiXi;S2-X2\SX){eO;SX\JX)D'};^j{*{p) , 

where D{p) is a Wigner-Z? matrix and 5'i_2 and Ai^2 are 
respectively the spins and hehcities of hadron 1,2. The 
spins are coupled to S = Si (>^ S2, ^ is the partial wave, 
J ~ £ (E) S is the total angular momentum and M is 
its z component. However, in all but the simplest cases, 
multi-hadron operators constructed by subducing Eq. Q 
into irreducible representations of the lattice symmetry 
can mix single-hadron operators transforming in differ- 
ent lattice irreps. Therefore, we prefer Eq. ([3]) where 
such mixings do not occur. The single-hadron operators 
transforming in definite lattice irreps can be optimised 
variationally, as shown for the pion in Section |IV[ 

Here we concentrate on the operators to be used to 
study two-pion states; the generalisation to other multi- 
hadron states is given in Appendix [A] The flavor struc- 
ture of the operators, for example the projection of tttt 
onto definite overall isospin, /, determines which combi- 
nations of Wick contractions appear in the calculation of 



the correlators (Section IH I . Because this flavor structure 



generally factorises from the spin and spatial structure we 
will not discuss it in detail here. However, because we are 
considering two identical pions, Bose symmetry requires 
the overall wavefunction to be symmetric under the in- 
terchange of the two pions. Therefore, in the 1 = 2 case 
we are considering here or / = 0, the symmetric flavor 
part requires a symmetric spatial part (even partial waves 
with positive parity). In contrast, 1 = 1 requires an an- 
tisymmetric spatial piece (odd partial waves with nega- 
tive parity). In addition, these operators have definite 
charge-conjugation parity, C = -1-1, for neutral combina- 
tions, generalising to G-parity for charged combinations; 
for brevity, in the following we omit the C-parity labels. 



A. Single-hadron operators 

Respecting the reduced symmetry of a finite cubic lat- 
tice, the J = 0~ pion at rest suhduces onto the one- 
dimensional A^ ~ A^ irrep of the double-cover octa- 
hedral group with parity, O^ . In Refs. [TTl [H] we dis- 
cussed how operators with a definite continuum J^ and 
Jz-component M, O'^ '^'^ {k = 0), can be constructed out 
of fermion bilinears featuring gauge-covariant derivatives 
and Dirac gamma matrices; the extension to baryons was 
described in Ref. [2]. The appropriate lattice operators 
were formed by subducing these continuum operators into 
octahedral group irreps. Table In] summarises how differ- 
ent integer continuum J subduce into octahedral group 
irreps - here we focus on the irreps relevant for mesons 
but the discussion applies equally to the irreps appro- 
priate for half-integer spin. In the case of a J^ = 0^ 
operator subducing to A^ = A^ this subduction is triv- 
ial, 

C)[o:l(0)= 00^(0) . 



p 


LG(P) 


A^ 


r 






A^ 


0±, 4±, ... 






Tt 


1±, 3±, 4±, ... 


[0,0,0] 


0? 


Tt 


2±, 3±, 4±, ... 






E^ 


2±, 4±, ... 






At 


3±, ... 



P 


LG(P) 


A 


|A|(''' 






Ai 


0+, 4, ... 






A2 


0-, 4, ... 


[0,0, n] 


Dic4 


E-2 


1, 3, ... 






Bi 


2, ... 






B2 


2, ... 






Ai 


0+, 2, 4, ... 


[0,n,n] 


Dic2 


A2 


0", 2, 4, ... 
1, 3, ... 






B2 


1, 3, ... 






^1 


0+, 3, ... 


[n, n, n] 


Dies 


A2 


0-, 3, ... 






E2 


1, 2, 4, ... 


[n, m,0] 


C4 


Ai 


0+, 1, 2, 3, 4, ... 


[n, n, m] 


A2 


0-, 1, 2, 3, 4, ... 



TABLE II. The pattern of subductions of the continuum spin, 
J < 4, (for P = 0) and heUcity, |A| < 4, (for P / 0) into 
lattice irreps, A l2^. Here f\ = P{—1)'\ P is given in units 
of ^ and n, m are non-zero integers with n ^ m. We show 
the double-cover groups but only give the irreps relevant for 
integer spin. 



At non-zero momentum, fc, the symmetry is reduced 
further: the relevant symmetry group is the little group, 
the subgroup of allowed transformations which leave k 
invariant [25] . In an infinite volume continuum the little 
group is the same for each fc; with only the constraints 
arising from rotational symmetry, states are now labelled 
by the magnitude of helicity, jAj, rather than J. On a fi- 
nite cubic lattice with periodic boundary conditions the 
allowed momenta are quantised, k = "^{n^m^p) where 
n,m,p are integers, and in general there are different 
little groups for different types of momentum. We de- 
note the little group for k by LG(fc) and for convenience 
define LG(()) = O;, . The pion subduces onto the one- 
dimensional A = A2 irrep of the appropriate little group 
(at least for all jA^j^ < 14 (^)^). Table |ll] shows the 
pattern of subductions of the helicities inro the little 
group irreps. In Ref. [27! we presented a method to 

construct subduced helicity operators, O)^ (k), and 

showed that these are useful for studying mesons with 
non-zero momentum on the lattice. For a J^ = 0^ oper- 
ator subduced into the A2 irrep the construction is again 
trivial, 

oi4°;'"i(fc) = o°"(fc). 

When we use the variational method to find the opti- 



mal linear combination of operators to interpolate a pion, 
we will include in the basis operators of other J subduced 
into A^ (for fc = 0) and other helicity A subduced into 
A2 (for fc 7^ 0). The pattern of subductions is given in 
Table In] the subduction coefficients for zero momentum 
are given in Ref. [12j and those for non-zero momentum 
are given in Appendix |A[ Henceforth, we will use 7r(fc) as 
a shorthand to represent O^- {k = 0) or O^j {k ^ 0) as 
appropriate. 



B. Multi-hadron operators 

In general, a tttt creation operator can be constructed 
from the product of two single-pion creation operators, 



, Jfcl,'i:2 t 
^ 'P,A,ii 



J2 C(P, A, m; fci;fc2)7rt(fci)^^(fc2), 



feiG{fei}* 

fc2G{fe2}* 

ki+k2=P 



(5) 
where 7r(fc) is a single-pion operator and C is a Clebsch- 



A- of 



Gordan coefficient for Ai (g) A2 — >■ A with Ai 2 
O^ if fci,2 = and Ai,2 = ^2 of LG(fci,2) if ^1,2 ^ 0, 
and where A is an irrep of LG(P). For present purposes, 
the particular construction of 7r(A:) from quark and gluon 
fields is not important. It is only necessary that 7r(/c) 
transforms in the appropriate lattice irrep. 

The sum over fci 2 is a sum over all momenta in the 
stars of ki_2, which we denote by {ki^2}* , and by which 
we mean all momenta related to fci^2 by an allowed lattice 
rotation. In other words, the sum is over i?fci 2 V i? G 
0/j ; the restriction that ki + k2 ^ P is equivalent to 
requiring R £ LG(P). We will write ki, ^2 and P in units 
of ^, using square braces to indicate the suppression 
of the dimensionful factor, i.e. P = [1,0,0] denotes a 
momentum of ^(1,0,0). 

The Clebsch- Gordan coefficients, C, can be determined 
by a group theoretic construction. When P = ki — 
^2 = 0, there is only one momentum direction in the 
sum and C arc just the usual Clebsch- Gordan coefficients 
for 0,j [28 . In the case of two pions the only relevant 
Clebsch-Gordan is the trivial A^ (g) A^ — > A^ , C = 1, 
giving a two-pion operator in the A'l irrep. 

For the two-pion system with fci ^ but overall at 
rest, P = ki + k2 = 0, k2 = — fci, the Clebsch- Gordans 
required are those for A2{{ki}*) (81 ^2(1^2}*) ~^ A^ with 
A2 of LG(fci) and A^ of O^ . The irreps, A ^, ar ising are 
given in Ref. [2S] and summarised in Table III We dis- 



cuss how to calculate the corresponding explicit Clebsch- 
Gordan coefficients using the induced representation and 
give values in Appendix [A} 

For the remaining case, P ^ 0, we require the Clebsch- 
Gordan coefficients for A2{{ki}*) ® ^2(1^2}*) — >■ A, or 



p 


fci 


fc2 


A(P) 






[0, 0, 0] 


[0,0,0] 


At 




[0,0,0] 

oP 


[0,0,1] 


[0,0,-1] 


A+,E+, (T-) 




0,1,1] 


[0,-1,-1] 


A+,T+,E+, {T-,T^- 


) 




1,1,1] 


[-1,-1,-1] 


A+,T+,{T,-,A^) 






0,0,2] 


[0,0,-2] 


At,E+, (TD 






0,0,0] 


[0,0,1] 


Ai 






0,-1,0] 


[0,1,1] 


Ai, E2, Bi 




[0,0,1] 
Dic4 


[-1,-1,0] 

[0,0,-1] 

[0,-1,-1] 


[1,1,1] 

[0,0,2] 
[0,1,2] 


Ai, E2, B2 

Ai 

Ai, E2, B\ 






-2,0,0] 


[2,0,1] 


^1, E2, Bi 






-1,-1,-1] 


[1,1,2] 


Ai, E2, B2 






0, 0, 0] 


[0,1,1] 


Ai 






[0,1,0] 


[0,0,1] 


Ai, (Bi) 






-1,0,0] 


[1,1,1] 


Ai, B2 




[0,1,1] 


[1,1,0] 


[-1,0,1] 


Ai, A2, (Bi, B2) 




DiC2 


0,1,-1] 


[0,0,2] 


Ai,Bi 






0,-1,0] 


[0,2,1] 


Ai, Bi 






[1,-1,1] 


[-1,2,0] 


Ai, A2, B\, B2 






[1,-1,0] 


[-1,2,1] 


Ai, A2, Bi, B2 






0,0,0] 


[1,1,1] 


Ai 




[1,1,1] 

DiC3 


1,0,0] 


[0,1,1] 


Ai, E2 




[2,0,0] 
[1,-1,0] 


[-1,1,1] 
[0,2,1] 


Ai, E2 
Ai, A2, 2E2 






[-1,0,0] 


[2,1,1] 


Ai, E2 





TABLE III. The two-pion operators for each P; also shown is 
LG(P) - we show the double-cover groups but only give the 
irreps relevant for integer spin. Example momenta fci and fc2 
are shown; all momenta in {fci}* and {fc2}* are summed over 
in Eq. Isl Swapping around fci and ^2 gives the same operators 
up to an overall phase. The irreps given in parentheses do 
not occur for two identical bosons with a symmetric flavour 
coupling (e.g. tttt in 7 = or 2) because of the constraints 
arising from Bose symmetry. 



if fc2 — (), A2{{ki}*) ® A^ ((}) — > A and correspondingly 
for fci — Q. Again, these are calculated using the induced 
representation as discussed in Appendix IX] and we give 
the irreps which arise in Table ]lTT] 



In this work we restrict ourselves to P = [0,0,0], 
[0,0,1], [0,1,1] and [1,1,1], and the vario_us combina- 
tions of fci and fc2 used are given in Table 



III 



Because 

the two pions are identical bosons, Bose symmetry re- 
quires them to be symmetric under interchange and we 
only use operators with the correct symmetry for isospin- 
2. However, for completeness, those operators with the 
wrong symmetry are shown in parentheses in the table. 

We want to use these operator constructions at both 
the source and the sink of correlation functions. This re- 
quires us to be able to project single-pion operators onto 
a given momentum at arbitrary times, something that 
can be achieved efficiently using the distillation method- 
ology |5D] . 



III. DISTILLATION AND CORRELATOR 
CONSTRUCTION 

Within distillation [3D] , we construct operators capable 
of interpolating a single pion of momentum k from the 
vacuum as 



TT\k,t) 



\ "^ ik-x 



i,u, 






^ 



{x,t), 



(6) 



where the Tt are, in general, operators acting in space, 
color and Dirac spin on a time slice, t, whose explicit con- 
struction is described in detail in Ref. ^7}. The quark 
fields tp in Equation [6] are acted upon by a distillation 
smearing operator Dg. that emphasizes the low momen- 
tum quark and gluon modes that dominate low mass 
hadrons. This smearing operator is defined as 



A^v, 



o'Ji^,y;t) 



e'^''"/'CMt)aHy,t) (7) 



where A„, ^^(a?, i) are the n*^ eigenvalue and eigenvector 
(in color, i, and position, x) of the gauge- covariant three- 
dimensional Laplacian operator on a time-slice, t. In the 
present study, the smearing weight a is set to and the 
number of vectors used is iVvocs = 64, 128, 162 on the 
L/fls = 16, 20, 24 lattices respectively (a shorthand D is 
used to represent Dj^^o)- 

The outer-product nature of the distillation smearing 
operator is such that correlators can be factorized into 
products of factors containing only propagation and fac- 
tors containing only operator construction. The propa- 
gation factors, T (called "perambulators"), and momen- 
tum projected operators, $, are constructed as matrices 
in the space of the eigenvectors (the distillation space) : 
where Tr„n{t',t) = ^l{t')M-^t' ,t)^^{t) and $„„(i) = 
^l^{t)Tt^rn{t), and M is the lattice representation of the 
Clover-Dirac operator for the light quarks used in this 
study. 

As outlined in Section |TTj two-hadron operators are 
constructed from sums over products of two single- 
hadron operators of definite momentum, as in Eq. [5] 
The resulting correlators for the tttt operators are of the 
generic form 



Q,(t',t) = (oiM^(i')-M;(t)|o), 



(8) 



where each operator tt is of the bilinear form given in 
Equation [61 For isospin-2, quark integration leads to 
only Wick contractions featuring quark propagation from 
source time t to sink time t'; there are no annihilation 
contributions. The resulting traces are over the set of 
eigenvectors used in Equation [7] which is much smaller 
than the full lattice space, allowing for the efficient com- 
putation of the correlation functions. In particular, it 
is the factorization of the smearing operator that allows 
for the projection of both the source and sink operators 
onto definite inter-pion momentum, something that is not 



possible in the traditional "point-all" method. This fac- 
torization allows for the construction of the full hermi- 
tian correlation matrix among source and sink operators 
in Eq. [8] and hence makes possible the application of 
the variational method }8t -tl0 . In this method, the mani- 
fest orthogonality among states provides the essential key 
for determining high lying excited states and separating 
nearly degenerate states. 

To increase statistics, the correlation functions in Eq.[8] 
are averaged over multiple time sources. The number of 
time sources, along with the number of eigenvectors of 
the Laplacian, N^ecs, and the number of configurations 
for each of the three volumes used in this study are shown 
in Table [H 



IV. OPTIMISED PION OPERATORS 

In our previous study of tttt isospin-2 scattering [7| we 
made use only of the simplest composite QCD operators 
capable of interpolating a pion, ^ -ipD^j^O^tp (Eq.p[with 
r = 75) where the distillation smearing operator Do- in 
Eq.[7]took on two different values of the smearing weight 
a. As well as interpolating the ground-state pion from 
the vacuum, this operator has significant amplitudes to 
interpolate various excited mesons with pion quantum 
numbers {tt*). In correlation functions, the contribution 
of the excited states will die away more rapidly than the 
ground-state (see the decomposition in Equation [ij , but 
at modest times, the excited states are present to some 
degree, as shown in Figure [l] For consideration of tttt 
scattering, these excited-state contributions are an un- 
wanted pollution in our correlators that ideally we would 
like to be absent. Their presence forces us to analyse 
TTTT correlators only at large times where effects of the 
finite-temporal extent of the lattice are more keenly felt. 

In principal if we could find an operator which has in- 
creased overlap onto the ground-state pion and reduced 
overlap onto low-lying excited states, its use would lead 
to TTTT correlators that are truly dominated by tttt at 
smaller times, with the contribution of unwanted tttt* 
being reduced. Our approach to finding such an "opti- 
mised" single-pion operator is to variationally diagonalise 
a matrix of single-hadron correlators in a basis of opera- 
tors, taking as our optimised operator the linear combi- 
nation of basis operators having lowest energy. 

The basis of operators used is as described in Section 



II A and presented in detail in Refs. [TT[ [T^ ^\ . It cor- 
responds to fermion bilinears with Dirac gamma matri- 
ces and gauge-covariant derivativeaj between them, con- 
structed to be of definite spin or helicity in a continuum 



* in this work we use all operators with the correct quantum num- 
bers constructed from any possible gamma matrix and up to 
three derivatives (for an operator at rest) or up to one derivative 
(for an operator at non-zero momentum) 



theory and then subduced into irreducible representa- 
tions of the octahedral group or the appropriate little 
group. For the pion this is A^ for zero momentum and 
A2 for all the non-zero momenta that we consider. 

The variational analysis corresponds to solution of the 
generalised eigenvalue problem IBHTU] 



C(t)«(") - Xnit)C{to)v^''^ 



(9) 



where the state energies are obtained from fits to An (t) ^ 
g-£;„(t-to)^ The optimal combination of operators, Oi, to 

interpolate state |n) from the vacuum is ilj, = X^i '^i" ^I- 
Our implementation of the variational method is de- 
scribed in Ref. [H]. 

In Figure [2] we show, for a range of momenta, the im- 
provement obtained using the "optimised" pion operators 
alongside the simple ipO^^Dip operators, where clearly 
the correlators computed with the optimised operators 
relax to the ground state more rapidly that the simpler 
operators, typically at or before lOaj from the source (a 
time comparable with the values of tg found to be optimal 
in solution of equation ^ . 

Use of these optimised operators will lead to some con- 
fidence when dealing with tttt correlators where for times 
> lOat away from the source, we will be able to largely 
neglect the contribution of tttt* states. 



V. PION MASS AND DISPERSION RELATION 

As well as the volume dependence of energies of multi- 
hadron states owing to hadron interactions suggested by 
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FIG. 1. Contributions of ground state (n — 0) pion (red) and 
excited pion states (other colors) to the single pion correlator 
at zero momentum, C{t) = ([ija-ysDil^) (t) ■ (^□75D7/))(0)) 
and A'vecs = 162 on the 24^ lattice. Summed contribution of 
all states indicated by the grey curve. Excited state pions are 
observed to contribute significantly until t > 20 at . (Excited 
state contributions determined from the results of variational 
analysis using a large operator basis, see the text) 
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FIG. 2. Effective masseajof single-pion correlators computed 
using ipD'y^a^l; (darker shades, squares) and "optimised" op- 
erators, f2n=o (lighter shades, circles). Shown for a range of 
momenta on the L/as = 24 lattice. 

^ throughout this paper we define the efi'ective mass of a 

Zat 



correlator C(i) to be rrtcff 



yC(t + Zat) 



the Liischer formalism, there can also be exponential de- 
pendence of single- hadron energies on L. We can attempt 
to determine any such behavior for the pion by computing 
its mass on the three volumes at our disposal. In Figure 
[3] we show the pion mass extracted on our three lattice 
volumes where there is seen to be very little volume de- 

pendence ( ZltL/Zl\% -0-990(4)). In 1^, NPLQCD 
suggest a xPT motivated form for the L dependence, 



^ — rriTrL 



m-^[L) = 7n„ + c 



(rri^L) 



3/2" 



(10) 



Fitting to this form we find atrriT^ = 0.06906(13) and 
atc = 0.24(10) in good agreement with NPLQCD's 
0.069073(63) (62), 0.23(12)(7) respectively. We use 
OtTTiTr = 0.06906(13) as our best estimate for the pion 
mass in all subsequent calculationf[j 

A complication which arises from our use of an 
anisotropic lattice is the need to determine the precise 
value of the anisotropy, ^, which relates the spatial lat- 
tice spacing as to the temporal lattice spacing at = fts/C- 
The anisotropy appears in the dispersion relation of a 
free-particle, where the periodic boundary conditions in 



space lead to allowed momenta p ~ 



27T 



integer n^ 



I'y, 



SO that 



'^x^ ^y 5 ^2 



{atEn^f = {atm) + — 



2 1 / 27r 



e \L/o 



for 



(11) 



^ fitting the same data to a constant leads to atrriT^ = 0.06928(18) 
with a somewhat poorer fit, X^ l^dot = 3-0- 



if we assume that mesons on the lattice have a 
continuum-Hke dispersion relation. Whether this is a 
good description will be determined by explicit fits to 
extracted pion energies at various momenta. In Figure HI 
we show pion energies on the three volumes along with 
a fit to Equation [TT] with ^ as a free parameter. The 
fit is acceptable leading to ^ = 3.444(6). Using other 
parameterisations of the dispersion relation (adding a p* 
term, using cosh/sinh etc.), lead to fits which are indis- 
tinguishable within the thickness of the line in Figure |4] 
and to compatible values of ^. In the remainder of the 
paper we use S, = 3.444(6) as our best estimat^ 



VI. EFFECTS OF FINITE TEMPORAL 
EXTENT 

Our extractions of finite- volume tttt energy spectra fol- 
low from analysis of the time-dependence of correlation 
functions and the form of these time-dependencies is af- 
fected by the finite temporal extent of the lattice. The 
size of finite-T effects are generically determined by the 
size of e"™""^, which while small on these lattices, is large 
enough for its effects to be visible, particularly in the tttt 
sector. 

As an explicit example of a systematic effect whose 
origin will turn out to be the finite temporal extent 
of the lattice, we show in Figure [5] the effective mass 
of a very simple "tttt" correlator. The same "tttt" op- 
erator, J23 [V'^TsOi/'] (x) ■ Y.y ['/'°75ni/'] (y) appears at 
source (igrc — 0) and sink (t) . The effective mass of the 
raw correlator is observed to continue falling after ap- 
pearing to briefly plateau near an energy equal to twice 
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FIG. 3. Pion mass as a function of lattice spatial volume. 
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the pion mass. This behavior can occur if the correla- 
tor features, as well as a sum of exponentially decay- 
ing time-dependencies corresponding to discrete energy 
eigenstates, as in Equation [ll also a contribution that is 
constant in time. Such a term can be eliminated by con- 
sidering the shifted correlator, Cst{t) = C{t) — C{t -\- St). 
An effective mass of this construction with St = Sat is 
also shown in Figure [S] where it is observed to plateau 
to an energy slightly above twice the pion massF] A di- 
rect estimate of the size of the constant term comes from 
fitting C{t) to the form 



E 



AnC 



'E„t 



(12) 



where {An}, {^^n} and c are the fit parameters. A fit 
to the raw correlator over a time region 15 — ?> 43 with 
two exponentials and a constant gives a x^/^doi = 0.7 
and atEo = 0.13966(4), Aq = 8.76(8) x 10* and c = 
26.4(13), indicating a statistically significant constant 
term. We propose that the origin of the constant term is 
the finite temporal extent of the lattice and notice that 
2^06""'"^ = 25.4(2) is very close to the fitted value of 



® in correlated fitting to obtain atm-n and ^ simultaneously we 
find a relatively small correlation between the parameters and 
for error propagations in the remainder of the calculation we 
treat them as independent variables. 



"^ the very slow relaxation to the plateau is mainly due to not using 
optimised pion operators in this construction. 
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FIG. 5. EfTective masses of a "tttt" correlator as described in 
the text. Raw correlator (red squares) and shifted correlator 
(green diamonds). 



c. In the remainder of this section wc will attempt to 
describe the effect of finite-T on om" computed tttt corre- 
lators at rest and in-flight. 



A. Finite-T effects for correlators at rest 

Let us begin by considering a correlator constructed 
using pion interpolating fields of definite momentum, 

(jk[,k^,. = Uz (t)nZ (t) ■ nt (0)7r+ (0)), 

kiM ^ ' \ k{^ ' k'^^ ' ki^ ' k.2^ '/' 

where in this section the operator 7r+ interpolates a posi- 
tively charged pion from the vacuum. In practice we will 
always project these products into definite little group 
irreps, A, for a given P = fci -I- fc2 = ^i + ^2 ^s described 
in Section [nj With anti-periodic boundary conditions in 
the finite time direction, two-point correlators have the 
decompositiorjj 



C{t) ^ {0'{t)0{Q)) 

= tr[e-^^O'(i)O^(0)]/tr[e-^^] 



E< 



'EnT {En-E,.a)t 



(n|C(0)|m)(m|O^(0)|n), 



(13) 



in terms of eigenstates of the Hamiltonian, H\n) ~ En\n), 
which will be discrete in a finite spatial volume. The 
contribution to this sum we are interested in is the only 



* see 1311 for a discussion of these finite-T effects on the spectrum 
of single particle systems and [32J for discussion of many-hadron 
states. 



one to survive in the limit T -^ oo and is of the form 

'^(0|7r-7r-|(7r+7r+)n)((7r+7r+)n|7r+7r+|0)e-^"'* 



n 



E(^-)' 



-E^ 



'\U) 



where |(7r+7r+),-,) are 1 = 2 eigenstates. 

At finite T there are other terms in the sum in Equa- 
tion [13) the largest being of form 



E' 



-E^{p)T/-\- 



(vk,,w-£,wi4)(4l4(o)4^(o)K,-) 



5ri r Sr, r 

k^,k2 fc2,fciJ 



which has a time-dependence of 

z-r Zt- \dr, r Sri r 

ki k2 L k[,ki k'^,k2 

X \^'E^(k[)T ^~{E^(k'^)^E„(k\))t 

_^g-B,(fe^)Tg-(B,(fei)--Ex(fe^))tl (15) 

where zr = (vr- tt- 0). As a first example, consider the 

case of correlators in the tttt rest frame, P = 0, C^'^ ^{t)^ 
where this term becomes 



2 {^k) 



4 g-EAk)T 



(16) 



which is simply a constant in time. 

We may now address the observation made at the 
start of this section that the correlator constructed with 
k = k' = [0, 0, 0] has a clear constant term. Our 
analysis above suggests that its magnitude would be 
2 (z[ooo]) e"™""^, while the leading T-independent term 

is of form IZttJ) e~^^^*. In the limit of weakly- 
interacting pious we would expect Zl^J — > (z[ooo]) and 
as such c — >■ 2AQe~'^^^ . This appears to hold true to a 
rather good approximation in the data. 

We expect other finite-T terms to be negligibly small 
in practice; in particular a term often considered in 
single-particle analysis, Z'^e~^'^e^* , which turns expo- 
nentially decaying time-dependence into cosh-like time- 
dependence can be ignored here. It is suppressed by at 
least e"^™'"'^ and only becomes relevant close to i = T/2 
while we consider correlators only at earlier time^ 

In Figure [6] we show further evidence for the presence 
of the constant term in tttt correlators. These correla- 
tors, evaluated on the L/as = 24 lattice, using opti- 
mised pion operators, have P = [0, 0, 0] and use ttj: tt^ 



® In practical terms, while the constant term could contribute ~ 
10% of the correlator at t = 48, the extra term "in the cosh" 
would only be at the 1% level. Other contributions featuring 
(7r+7r+7r^|7r+7r+|7r='=} formally appear at 0{e~"^^'^), but their 
i-dependence ensures that they provide negligible contributions 
to the correlators. 
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FIG. 6. Fits to diagonal tttt correlators with P = [0, 0, 0] 
using the lowest allowed \k\ that gives rise to irrep A^. Cor- 
relator is plotted via e^^"'-'''* C(i) such that in the limit of 
non-interacting pions and T — ^ oo we would have a horizon- 
tal line. The solid red line shows the result of the fit using 
equation (121 while the orange dashed line shows the result of 



excluding the constant contribution, which should correspond 
to the T — > oo behavior. Fit parameters given in Table |IV[ 



products projected into definite irreps A-'^ constructed 
from the lowest allowed |fci|, \k2\ as detailed in Section 
Il](S+ ^ fc= [1,0,0], T+ ^ k = [1,1,0]). The results 



of the correlator fits (of the form given in equation ( 12 )) 
are presented in Table |IV[ where we see that the size 
of the constant term is in rather good agreement with 
4Aoe~^'''-'^^"^, the value in a non-interacting theory (see 
the Clebsch-Gordan coefficients in Appendix |A] for the 
appropriate combination of pion momenta) . Clearly the 
"polluting" constant term plays a significant role in the 
correlator as early as t ~ 25at and if we want to use 
timeslices beyond this point in variational analysis, we 
will need to take some account of its presence. 

Our solution is to completely remove the effect of all 
time-independent terms from correlators, by instead of 
considering C{t), using shifted correlators. 



Cstit)^C{t)-Cit + 6t), 



(17) 



which exactly cancel contributions constant in time for 
any choice of St ^ 0. The desired tttt contributions, equa- 



irrep. k x'/Ndoi Eo/2E^{k) 



~E^ 



[1,0,0] 
[1,1,0] 



0.9 
0.9 



10 



1.0014(17) 7.9(5) X 10" 
1.0002(16) 6.9(11) X 10" 



7.8 X 10" 

4.9 X 10" 



TABLE IV. Fits, using two exponentials in equation (12 1, 
to diagonal tttt correlators with P = [0, 0, 0] using the lowest 
allowed |fc| that gives rise to irrep A^. 



tion ( 14 ) , are changed only by a rescaling of the Z"^ to 

1/2 



TTTT TVTT 



1-e 



-El, St 



(18) 



This is just a change in scale of overlaps that for a given 
state, n, is common to all operators. Shifting then does 
not violate any of the conditions for carrying out a vari- 
ational analysis and we can proceed with use of Cst (t) in 
equation ([9]) to yield the finite-volume energy spectrum 



E^-R- 



B. Finite-T effects for correlators in-flight 

The unwanted contributions to correlators "in-fiight" 
(P 7^ 0) are not time-independent and cannot be removed 
by simply shifting in time. Following equation (13 1, they 
take the generic form 

X [^-EAk[)T^-{E^{-k[+P)-E^{k[))t 

_^^-E^i-k[+P)T^-{E^{k[)~E^{~k[+P))tl 

where the contributions of largest magnitude occur if ei- 
ther k[ or —k[+ P are equal to zero as then the finite-T 
suppression factor is only e"™""^. The largest "polluting" 
term in this case would not be a constant but rather have 



a time dependence 



-AE^t 



where AEj^ is the energy 



gap between a single pion of momentum k and one with 
momentum P — k. In the case P = this reverts to a 
constant in time as expected. 

Consider the concrete example of a correlator with 
7ri(,g|7r^jjp, at the source and 7rrQQg|7r,~pg, at the sink. In 
this case, as well as the desired term which is approxi- 
matel\F°l 



(.[ooo,)^(^[roo,)^"("-^^i^™'^ 



(19) 



we would have "polluting" terms 



(^[OOO]) (^[100]) {<^ 



-7n^T-{E^ 



100] _ 



_^ ^__BU001-r^_(„^_£;[1001)j 



(20) 



^'^ this would be exact for non-interacting pions, in tttt 1 = 2 scat- 
tering the interaction is weak so the approximation should be a 
reasonable guide. 
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FIG. 7. Simulated contributions to a correlator (fci = 
[0,0,0], fc2 = [1,0,0]) of the desired (T ->• oo, red) term 
polluting" (finite-T) terms from equa- 



equation 
tion 



20 



19 



and two 
the first term (leading, green dashed) and the sum 
of the two terms (leading plus subleading, blue, dot-dashed). 
Observe that in the time region we will consider, the leading 
term dominates over the subleading term. 



where as shown in Figure [7j the first of these pollut- 
ing terms is expected to dominate the pollution for the 
time regions we consider. We can observe the effect of 
this leading pollution term in fits to correlators having 
P = [1, 0, 0] computed on the L/ag = 24 lattice using op- 
timised pion operators - in Figure [8] we show the irreps 
A = Ai, Bi, B2, constructed using the smallest allowed 
magnitudes of pion momentum. The fit form (which ne- 
glects the subleading pollution) is 



A, 



-E„t 



ce 



^{E^(k„ 



.)-E^ik^,„))t 



(21) 



with {^n}, {^-n} and c as fit variables, using fixed ET^{k) 
obtained from the dispersion relation (equation 11). 



It would appear that these diagonal correlators can 
be reasonably well described by the fit form proposed 
indicating a small but statistically significant impact of 
finite-T effects on the correlators. We will need to address 
these terms in any variational extraction of the in-flight 
TTTT spectrum. Our approach is to remove the worst of 
the pollution exactly and settle for approximate reduc- 
tion of less acute terms. The largest polluting term has 



a time-dependence oc e 



-E„(k^ 



.)Tg-AB„ 



where k,„;„ is 



the lowest momentum that appears in any of the cor- 
relators making up our correlator matrix, and Aii^min is 
whatever positive energy gap appears in the correspond- 
ing time-dependence. This term can be converted into a 
constant by forming the following weighted correlator. 



C{t) 



,AB,„i„ t 



C{t), 



(22) 



and the constant term can be removed by then shifting 
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FIG. 8. Fits to diagonal mr correlators with P — [0, 0, 1] 
using the lowest allowed |fci|, |fc2| that gives rise to irrep A. 
Correlator is plotted via e^^^'^i'+^'f^^))* ^(t) such that in 
the limit of non-interacting pious and T — >■ 00 we would have 
a horizontal line. The solid red line shows the result of the fit 



using equation (211 while the orange dashed line shows the 



result of excluding the contribution proportional to c, which 
should correspond to the T — >■ oo behavior. 



the weighted correlator, 

Cst{t)^C{t)-C{t + St). 



(23) 



We refer to these as weighted-shifted correlators. The 
exact same weighting and shifting procedure is applied 
to every element of the matrix of correlators such that 
the effect of the weighting is to shift all energies down by 



a common AE'min- This can be corrected for by adding 
Ai?min to the variationally obtained spectrum. 

In summary, while finite-T effects are modest in our 
two-pion correlators, precision extraction of a tttt energy 
spectrum requires that we account for them in our analy- 
sis. Through appropriate weighting and shifting of corre- 
lators before applying the variational method, we beheve 
that we are able to remove the leading systematic effects 
leaving only sub-leading effects that we find to be smaller 
than our level of statistical uncertainty. 



VII. FINITE- VOLUME SPECTRUM 

We compute correlator matrices in each irrep P, A 
using the basis of operators defined in Section [llj Af- 
ter modifying the correlator matrix with the appropriate 
weighting and/or shifting as described in the previous 
section, the spectrum is obtained by solution of the gen- 
eralised eigenvalue problem, equation |9] Each irrep is 
considered independently and the entire procedure is re- 
peated on each of the three lattice volumes. The two-pion 
operators used are given in Table |V] and the numbe r of 
operators for each P and irrep are given in Table 



VI 



We illustrate the method here with the example of the 
P = [1,0,0], A = Ai irrep on the i/a., = 24 lattice. 



A. Example of P = [0, 0, 1], A = Ai 

As an explicit example of our variational fitting pro- 
cedure consider the P — [0,0, 1], Ai irrep evaluated on 
the L/tts = 24 lattice. Our basis of operators here is 
obtained by applying equation [s] and is thus of the form 



12 



tion expected to be present for t <to (our reported spec- 
tra are just the values of E, E' is discarded). The fits 
are very good and the absence of any significant upward 
curvature at larger t (as in Figure l7| suggests that our 
weighting-shifting procedure has removed the bulk of the 
finite-T pollutioi|3 

The solution of equation [9] also provides eigen- 
vectors v^"' which can be converted into overlaps. 



,[fcl,fc2]t 

'[001],Ai 



fei6{fci}* 

fc26{fe2}* 

fel-|-fc2 = [001] 



with constructions using the pion momenta given in 
Table Ivj This gives a correlation matrix, C{t), of di- 
mension 7. In order to remove the largest finite-T ef- 
fects, as discussed in the previous section, the weighted- 
shifted correlation matrix, C{t) is formed, using AEj^i^ — 
ii^^([0,0, 1]) — 77i7r. This matrix is analysed using equa- 
tion [9]- for to — 14at, the obtained principal correlators. 
An (t) are shown in Figure M along with fits of the form, 

-E{t-to) , A„-E'{t-to) 



X{t) = (1 - A)e- 



Ae- 



(25) 



where E, E' > E and ^ <C 1 are the fit parameters. The 
second exponential allows for the excited statq^ poUu- 



z^':K 



C{[0,0,l],A,;h;k2)TTHM)7T\k2) 



(24) 



Vt[fcl,fc2] 

,Ai 



[001], A Jo) '^'^i^g 
Our method of so- 



[fe^,fc,] - ((7r7r)„;[0,0,l],Ai|(7r^ 

[ki,k2] ^ ^ ' ' Iki,'=2] 

lution of the generalised eigenvalue problem treats each 
timeslice independently such that we actually obtain 
w("'(t) and thus Z{€). This time-dependence is fitted to 
a constant (or a constant plus an exponential if that is 
required to get a good fit) and the resulting constant is 
rescaled to undo the effect of the shifting of the correla- 
tors in the manner prescribed by equation |18[ 

The overall quality of description of the correlators by 



the variational solution can be seen in Figure 10 along 
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TABLE v. The two-pion operators used presented for each P 
on various volumes; also shown is LG(P). We give only the 
irreps that we considered in this work. Example momenta 
fci and fc2 are shown; all momenta in {fci}* and {^2}* are 
summed over in Eq. [5] 



^^ by "excited states" here we might have several types, including 
TTTT with large relative momenta, tttt* and other inelastic contri- 
butions. 



^^ such upward curvature is seen in variational analysis of the raw 
correlator matrix, C{t). 
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FIG. 9. Principal correlators from solution of equation lol applied to the weighted-shifted correlator matrix C{t) for P — [0, , 1] 



25 



A — Al -with to ~ 14at. Plotted is e^'*"*'"A(i) against t/at along with fits to the time-dependence according to equation 
Also plotted in the bottom-right are the effective masses of the principal correlators (with the energy weighting A^min corrected) 
and the fit values E superimposed as horizontal bands. All energies are those in the frame in which the lattice is at rest 
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FIG. 10. Diagonal elements of the weighted-shifted correlation matrix for P = [0,0,1], A = Ai: C[fej_g2](i) and their 

Plotted is e ^'" C{t). 
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TABLE VI. The number of two-pion operators used for each 
P and irrep on the various lattice volumes. 
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It is apparent in Figure 10 that the basis of operators, 
defined by equation [24] is rather close to a diagonalis- 
ing basis and this can be clearly seen in Figure [TT] which 
shows the Z values for each state n and each operator 
[ki,k2\- This indicates that the finite- volume tttt eigen- 
states are close to being states of definite pion momentum 
which agrees with the expectation that the 1 = 2 inter- 
pion interaction strength is weak and the observation of 
only small shifts from non-interacting tttt energies. It 
is interesting to note that the largest deviations from 
diagonal behaviour, i.e. the largest mixing of the non- 
interacting state basis, occurs for levels which are very 
close in energy. This is precisely what we would expect 
from perturbation theory, where small energy denomi- 
nators enhance mixing of near-degenerate states. That 
we are able to resolve this mixing with a high degree 
of confidence is an advantage of our use of a variational 
approach. 



B. Volume dependence of tttt spectra 

We perform this analysis procedure independently for 
each P, A on each volume. The energies obtained are in 
the frame in which the lattice is at rest, and can be more 
usefully expressed in the tttt center-of-momentum frame. 



Er, 



EL 



\P\- 



{citEa 



(flt-Eiat) ~ 52 ( 



27r 






1/2 



(26) 



where we use the anisotropy, ^, determined from the pion 
dispersion relation in Section |V] In Figures |12|13|14|15| 
we show the volume dependence of the extracted center- 
of-momentum frame energy spectrum along with the en- 
ergies of pairs of non-interacting pions carrying various 
allowed lattice momenta. 

In all cases we observe small energy shifts, with 
the largest shifts in Ai irreps, reflecting the expected 
strongest interaction in S'-wave scattering. 



all states for a given operator [kx , ^2] is unity. 



with an indication of how much each (7r7r)n) state con- 
tributes to each of the diagonal correlators. These con- 
tributions are reconstructed from the results of the varia- 
tional analysis by building the sum in equation [Ml state- 
by-state. The description, as one would expect, is excel- 
lent for t > to; indeed the ability to get a good description 
of the correlators using only the number of states equal 
to the basis size is our condition to determine an appro- 
priate value of to |33j . That we are able to countenance 
a value as low as to = lAat is due to our use of opti- 
mised pion operators so that tttt* contributions to the 
correlators are much reduced. 
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FIG. 12. Extracted center-of-momentum frame energy 

spectra for P — [0,0,0] irreps A^ , E'^ jT^ ■ Also shown (in 



red) are non- interacting pion pair energies, \/m'^ + \ki\'-^ + 



/ 



w.| + 1^2!^ whose uncertainty is determined by the uncer- 
tainty on atTTiTT and ^ determined in Section IV] Grey area 
represents opening of inelastic {4n) threshold. 
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FIG. 13. As Figure 12 for P = [0, 0, 1] 
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VIII. PHASE-SHIFTS FROM FINITE- VOLUME 
SPECTRA 

The formalism to relate the amplitude for two-particle 
elastic scattering in partial waves labelled by angular mo- 
mentum £, to the spectrum of states in a finite cubic spa- 
tial volume, is laid down in [2], with extensions to the 
case of a moving frame presented in [16til8j . Because 
we are considering tttt scattering in isospin-2 where only 
even (. occur, there is a one-to-one correspondence be- 
tween the irreps of the symmetry group relevant for the 
Liischer formalism in a moving frame and the little group 
irrepa^ and so we will refer to the little group irreps. 

The formalism can be compactly expressed in a single 
equation. 



det 



E(pc.)-U(^^^) 



27r / 



= 0. 



(27) 



U is a formally infinite-dimensional matrix of known 
functions whose rows and columns are each labelled by 
the pair (^, n), C/f„;f'n'- {^} are the angular momenta 
which subduce into the irrep. A, and n is an index in- 
dicating the n"^ embedding of that I into this irrep; the 
pattern of these subductions is given in Table |VII[ U is 

a function of the dimcnsionless variable q^ — (^^f^) , 
featuring the center-of-momentum frame scattering mo- 
mentum and the spatial length of the cubic lattice, L. 

E is a diagonal matrix, independent of L, which 
encodes the scattering amplitude through the elas- 
tic scattering phase-shifts, 5i,(pan)^ as Ein-Cn' = 

U is conveniently expressed in terms of a matrix M 
as U = {M + il){M-il)~'^ 
elements of M using 



where we can obtain the 



■^£;^;'^^(9')^A,A"5.,.' 



E 

A=±|A| 



E 

' = ±|A'| 



-.?y,A* 



iW 



i(p) 



sr:.. D'ZiR) ■ M)^U'^'i^) ■ sr. ,„ D^A-m 



f,,A' 



^in 



P,A,/i " 



P.A,^' ™'->^' 



(28) 



r 



In this equation, i? is a rotation carrying the J^ quan- we may be justified in making a finite truncation in L 

tisation axis (0,0, P) into P, with dU^{R) relating J^ 

values, m, to helicities, A. A convention for constructing 
R is given in 1271 . 5-' is the subduction from helic- 

ity A to the fi^^ row of the lattice irrep A (see Appendix 
[A| ). Different magnitudes of helicity, |A|, |A'| give rise to 
the different embeddings n,n'. The "reflection parity", 
i] = P{—lY = + for a system of two pseudoscalars. 

■^em-i'm' i^ ^^^ same object defined in equation (89) of 
PJ3j where it is expressed in terms of a known linear com- 
bination of generalised zeta functions of argument q^. 



One potential use of equation 27 is to take a scatter- 



ing problem where the amplitudes are known and find 
the corresponding spectrum of states in a certain finite- 
volume box. For a known set of scattering phase-shifts, 
{Si{pcm)}, the finite-volume spectrum on an L x L x L 
spatial lattice can be obtained by solving equation 27 
for discrete values of Pcm which give discrete values of 
-Eiat. Of course in practice, for any given lattice ir- 
rep, A, we need to truncate the infinite {£, n) basis to 
the set of phase-shifts {^^(pcm)} known to us. Fortu- 
nately, at low scattering momentum there is a hierarchy 
in di {pcm ) which follows from angular momentum conser- 
vation, (5f(pcm) ^ Pcm^^i such that Sq ^ 62 ^ 64 . . ., and 
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^'^ We note that the symmetry group relevant for the Liischer for- 
mahsm here is the subgroup of OjJ' under which P — > ibP rather 
than the constraint for httle groups that P —^ P. The irreps 
are similar to those of the little groups but have an additional 
"parity" label. 



TABLE VIL The pattern of subductions of / = 2 tttt partial 
waves, ^ < 4, into lattice irreps, A, where A^ is the number of 
embeddings of this £ in this irrep. This table is derived from 
Table Inl by considering the subductions of the £ for P — and 
the various helicity components for each £ for P 7^ 0. Here 
P is given in units of ^ and n, m are non-zero integers with 
n ^ m. We show the double-cover groups but only give the 
irreps relevant for integer spin. 



19 



A. A toy model of tttt scattering 

In order to demonstrate the formalism, we will briefly 
break away from analysis of lattice QCD obtained finite- 
volume spectra to consider a simple toy-model of tttt scat- 
tering in which the scattering amplitudes are known to 
us. The toy model is built from an effective range pa- 
rameterisation of elastic scattering in I = 0,2,4 partial 
waves. We have 






C0t(5(pcm) 



1 



1 



rtPlmi 



(29) 



2.5GeV^\ 



ai 2 

with parameters 

ao =:-0.8GeV"\ r^ 

a2 = -2.4GeV-^r2 = 0, 

aA = -5.0GeV~^ r4 eee 0, 

which happens to reasonably describe the experimental 
TTTT 1 = 2 scattering data up to a momentum p^m ^ 
0.6 GeV HH. 

Given this parameterisation and the choice 771,^ = 
0.396 GeV we solve equation [27] for the finite-volume 
spectrum in several irreps, (P, A), over a range of vol- 
umes, L = 2.0 -^ 5.0 fm. In Figure [T6| we show the center- 
of-momentum frame finite-volume energy spectrum for 
one example irrep P = [0, 0, 1], A = ^1. At each vol- 
ume we show the spectrum obtained from three different 
scattering parameterisations: the green squares show the 
spectrum with only 5- wave scattering {62 = S4 = 0), 
the blue circles include also D-wave scattering (^4 = 0), 
and the black diamonds correspond to all of ^0.2,4 being 
described by the effective range parameterisations given 
above. We observe that the contribution of higher par- 
tial waves to determining the finite- volume energy varies 
with excitation level. 

The problem to be solved in lattice QCD is actually 
the inverse of that just described - we start with the 
finite-volume spectrum determined through analysis of 
correlation functions and want to find the phase-shifts 
as a function of scattering momentum. If a given irrep 
received contributions from only a single i this would be 
relatively simple - we would solve equation [27] for un- 
known Sg{pcm) by inputting the determined value of Pcm 
extracted from E\3t using equation [26] The toy model 
construction indicates the potential difficulty with such 
a naive approach - equation [27] depends on the value of 
many Si simultaneously and on the face of it this is an 
underconstrained problem. 

Within the toy model we can explore the effect of the 
simplest possible assumption that higher partial waves 
contribute only negligibly - consider the spectr um in 
P ^ [0,0,1], A = Ai ior L = 3.5 fm. In Figure llTFa) 



we show the extracted Sq for the lowest four energy lev- 
els as a function of a supplied valu^^ of 62 (and with 



^4 = 0) . The naive assumption oi 62 = is seen to give 
reasonable estimates of (5o for the lowest two levels, but to 
be significantly discrepant for the next two levels. Vary- 
ing 62 between ±2|(52'^^'^'| (which we know because this is 
a toy model) gives the curves shown. Figure [l7) (b) shows 
the sensitivity to 64 assuming that 62 is known exactly. 

In this exercise we explicitly see that the influence of 
higher partial waves can vary significantly between levels; 
for n = 0,1,3 the infiuence of 62,4 is modest and given 
a "reasonable" estimate of their magnitude we could as- 
sign a systematic error on Sq that would encompass the 
exact result. On the other hand, no information can be 
obtained from level n = 2 without very precise knowledge 
of both 62 and ^4 at the corresponding scattering momen- 
tum. This will not be possible in any practical calculation 
and we must be careful to identify those cases where an 
energy level shows such extreme sensitivity. 

Thus even if our main aim is only to determine (5o(Pcm) 
we see that it is incumbent upon us to also estimate 
(5f>o. The easiest way to do this is to analyse the finite- 
volume spectra of irreps which receive no contribution 
from £ = 0, see Table |VII| Typically any irrep that fea- 
tures £ = 2 will also feature ^ = 4 so we have a similar 
problem of estimating S2 given no knowledge of S4 . For- 
tunately in the case under consideration where the inter- 
actions are weak we encounter situations in which energy 
levels in two different irreps have very similar energy val- 
ues. For example with P = [0,0, 1], the lowest level in 
E2 and the lowest level in Bi are both very close to the 
non-interacting (fci = [0,-1,0], ^2 = [0,1,1]) level and 
correspond to Pcm values of 0.03934, 0.03950 GeV respec- 
tively. In this case, to the extent that 52,4(Pcm) do not 
change significantly over the small difference in pcm , and 
the functions in U are not rapidly varying over the corre- 
sponding range in q^ , we can solve the coupled system of 
two equations 27 (for E2 and Bi) for the two unknowns 
62,64. This is demonstrated in Figure 18 |a) where the 
simultaneous solution of the two equations is seen to be 
reasonably close to the exact values. A similar extraction 



for two levels in E2, B2 is shown in Figure 18 ^b) 
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en- 



included in equation 27 in E as a fixed parameter 



Several level pairings of this type can be identified and 
an estimate of a few discrete values of ^2 7 64 can be made 
as shown by the purple points in Figure [T9] Fitting these 
points with an effective-range parameterisation, or using 
some other method to interpolate between the discrete 
points, we obtain our desired estimates of ^2j'54 for use 
in determination of Sq. 

The extracted values of 5q shown in Figure 
correspond to solving equation 27 for each 
ergy level in the Ai representation for P = 
[0,0,0], [0,0,1], [0,1,1], [1,1,1] with 62,64 fixed at our 
best estimate from interpolation between the determined 
62,64 points. The error bars indicate the uncertainty in 
60 obtained by varying 62,64 within an assumed 100% 
uncertainty. We observe that following such a procedure 
leads to a reasonable reproduction of the originally input 
toy-model phase-shifts. Note that we used only a single 
volume to obtain this result - using multiple volumes will 
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P=[0,0,l]L4i 



[0,-2,0], [0,2,1] 
[-1,-1,-1], [1,1,2] 

[0,-1,-1], [0,1,2] 

[-1,-1,0], [1,1,1] 
[0,0,-1], [0,0,2] 

-1,0,0], [1,0,1] 
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FIG. 16. Finite- volume spectrum for the toy model of effective-range parameterisations in the irrep P — [0,0,1], Ai. 
Green squares indicate the spectrum including only i = scattering, blue circles include £ = 0, 2 and black diamonds include 
i" = 0, 2, 4. Note that for many of the energy levels the squares, circles and diamonds he on top of each other. Red curves show 
non-interacting energies of pion pairs with momenta k\,k2. 
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FIG. 17. Lowest four energy levels (n = 0, 1, 2, 3 ) in toy model with volume L = 3.5 fm in irrep P — [0, 0, 1] , Ai. 

as a function of assumed values of S2 in range ±2|(52''^'^*| with 84 



(a) Sensitivity of 60 extracted from equation 

(b) Sensitivity of So extracted from equation 







as a function of assumed values of 84 in range ±2|(5|'"''''| with 62 = 62^'"'^ 



Boxes on far left indicate exact values of 60 at the corresponding scattering momenta. Arrows on x-axis indicate exact values 

of (52,4 



further improve the determination. 

An alternative approach to deaUng with the contribu- 
tion of higher partial waves is to parameterise all 61 {pcm ) 
one expects to contribute significantly in terms of a rela- 
tively small number of variable parameters. By perform- 



ing a global fit to all energy levels simultaneously, varying 
the parameters, one can attempt to find a description of 
the finite- volume spectrum that is best in a least-squares 
sense. Clearly in the case of this toy model, one could use 
the parameterisation given in equation [29] and by vary- 
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Pcm = 0.3941(8) GeV 

P= [0,0,1] 




Pcm = 0.5349(28) GeV 

P= [0,0,1] 




1.5 84/' 



FIG. 18. Simultaneous solution of two equations [27| for (52.4. Open black square shows exact values with uncertainty indicating 

the variation in (J™!''^' over the momentum region between the two determined pcm- 

(a) P = [0,0, 1], n = in E2 and n = in Bi. (b) P = [0,0, 1], n = 1 in E2 and n = in B2. 
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FIG. 19. Phase shifts, (So,2,4(Pcni) extracted from L = 3.5 fm spectrum using the method described in the text. Uncertainty in 
So indicates the effect of a conservative assumed uncertainty on ^2,4. Some points with very large uncertainty not shown. Toy 
model input phase shifts shown by the curves. 



ing parameters oq, tq, 02, 04 come to an exact description 
of the spectrum. We do not present this trivial result 
here, but we will return to this "global fitting" method 
in the next section when we consider the finite-volume 
spectrum obtained from lattice QCD computations. 
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B. Lattice QCD data 



We now return to consideration of the finite-volume 



spectrum presented in Section |VII| The first step in our 
"level- by-level" approach is to solve for ^2,4 using pairs of 
simultaneous equations [27] Pairs of levels below inelastic 
threshold that can be used to yield estimates for 62,4 are 



presented in Table VIII and are displayed by the filled 
points in Figures |20[ |21[ 64, is observed to be statisti- 
cally compatible with zero throughout the elastic region. 
There are also levels in irreps whose leading contribution 
is from i = 2 that do not pair and cannot be analysed 
using a simultaneous solution - these are considered in 
isolation, where the (small) role of ^4 is estimated and 
included as a systematic error, they are shown by the 
open points in Figure [20| 

Each of these S2 , S4 data sets can be described well by 
a scattering length fit, p^^^ cot Se{pcm) = 1/a^j and the 
resulting fit function is used to estimate the size of 62,4 at 
any pcm in the elastic region when determining Sq values 
from Ai irreps. As indicated in the previous subsection, a 
systematic error on So due to imperfect knowledge of 62^4 
is assigned by assuming a 100% error on the estimated 
values of 62^4. The resulting 60 points are displayed in 
Figure [22] where it is observed that the uncertainty from 
imperfect knowledge of 62,4 is typically much smaller that 
the statistical uncertainty. 
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FIG. 21. 5i values in elastic scattering region determined from 
finite-volume spectra. Filled points determined by simultane- 
ous solution of equations 27 (innermost errorbar statistical 
uncertainty, outermost errorbar reflects combined statistical 
uncertainty and uncertainty in atm^, ^ with errors added in 
quadrature) . 
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FIG. 22. So values in elastic scattering region determined 
from finite-volume spectra. Innermost errorbar is the statis- 
tical uncertainty, middle errorbar combined statistical uncer- 
tainty and uncertainty in (atm,r, C)i outermost errorbar re- 
flects total uncertainty including imperfect knowledge of ^2,4 
(all errors added in quadrature). Some points with very large 
uncertainty not shown. 



FIG. 20. ^2 values in elastic scattering region determined from 
finite-volume spectra. Filled points determined by simultane- 



ous solution of equations 27 (innermost errorbar statistical 
uncertainty, outermost errorbar reflects combined statistical 
uncertainty and uncertainty in atrriTy, ^ with all errors added 
in quadrature). Open points determined from single levels, 
with effect of ^4 estimated (innermost errorbar statistical un- 
certainty, outermost errorbar reflects combined statistical un- 
certainty and uncertainty in atm.^, ^ and 54). 



We now consider the second approach described above 
where the Se{pcm) are parameterised and by varying a 
small number of parameters a best description of all the 
finite volume spectra is obtained in a "global fit". Our 
procedure is to minimise a x^ with respect to the variable 
parameters in the parameterisation, which we denote col- 
lectively by {tti}. The x^ describes the similarity between 
the extracted finite-volume spectrum and the spectrum 
predicted by the parameterisation on the appropriate vol- 
umes, 
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L/as 



24 



levels 
[0,0,0], s+, 

[0,0,0], r2+, 



atPcm 



0.10766(23)(8) 
0.10764(23)(8) 



5-2/° 

-0.39(82)(67) 



-0.17(32)(22) 



24 



[0,0,1], Bi,n = 
[0,0,1], E2,n^0 



0.08427(25)(11) 
0.08418(25)(11) 



-0.40(47) (39) 



-0.05(26)(16) 
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[0,0,1], B2,n = 
[0,0,1], iJ2,n=l 



0.11412(29)(8) 
0.11393(28)(8) 



-1.60(80)(64) 



-0.78(69)(55) 



20 



[0,0,1], Bi,n = 
[0,0,1], E2,n^0 



0.10174(35)(9) 
0.10131(37)(9) 



-1.59(54)(36) 



-0.018(36)(17) 



TABLE VIII. Levels with very similar Pcm values used in simultaneous solution of equations |27| 
X\{a^}) - E E [p-(i; ^An) - pt\L; PAn; {a,})] C-\L; PAn; P'A'n') [p,^{L; P'A'n') - pt\L; P'A'n'; {aj) 



L PAn 
P'A'n' 



(30) 



Here we have p^^{L; PAn;{ai}) which is the particular 



solution of equation 27 which is nearest to Pcm(-^;^An) 
(with the parameters set to the particular values {oi}). 
The data covariance, C, accounts for the correlation be- 
tween determined energies computed on the same lattice 
configurations - different volumes correspond to indepen- 
dently generated lattice ensembles and hence are not cor- 
related. 

Statistical errors on the parameters, {ai}, are deter- 
mined by Ax^ — 1. Errors from the imperfect knowledge 
of atm^ and ^ are estimated by repeating the x^ min- 
imisation varying the mass and anisotropy within their 
respective errors. We treat these as independent system- 
atic errors, although they would naturally be reduced 
with increased numbers of gauge-field configurations at 
each lattice volume. 



"'1=0 

at.=2 



(- 



-4.13 ±0.07 ±0.06) • at 
-1.08 ±0.28 ±0.19) X 10^ 



1 0.5 
1 



xVA^dof = 121/47, 



where the quality of fit is insignificantly degraded. 
Clearly there is no need to invoke higher terms in the 
effective range expansion to describe the data. A fit to 
only those irreps where ^ = 2 is leading yields 0^=2 = 



(-1.51 ±0.31 ±0.21) 



10-^ 



with xVN^oi = 31/16, 



in reasonable agreement with the values obtained above. 
These various fits arc shown in Figure [23] along with the 
points determined using the "level-by-level" approach de- 
scribed earlier where good agreement between the two 
methods is observed. 



Fits with effective range and scattering length pa- 
rameterisations (equation 29 ) were attempted. These 



fits never indicated the need to include significant 
strength in the £ = 4 wave. A successful fit to 
all energy levels with an effective range parame- 
terisation oi i ~ and scattering length \\\ i — 2 
gives the following parameter values and correlations, 



ai^Q = (-4.45 ± 0.18 ± 0.06) • a* 
rfco =(-3.7 ±1.8 ±0.7) • a* 
afc2 = (-1.20 ± 0.29 ± 0.17) x 10^ 

xV^dof = 116/46, 



1 



0.9 0.4 
1 0.2 
1 



where the second set of uncertainties reflects variation 
of AjTOtt and ^ within their uncertainties. We see that 
the effective range in £ = is barely significant and 
is strongly correlated with the scattering length. The 
degree of correlation between i = Q and £ = 2 is mild. 
Given the lack of significance for rg, a fit with just a 
scattering length was attempted, yielding 
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FIG. 23. Upper: 5o{pcm) obtained through "global fits" 
to finite-volume spectra using effective range and scattering 
length parameterisations. Lower: 52(pcm) obtained through 
"global fits" to finite- volume spectra using a scattering length 
parameterisations. Also shown for comparison, the "level-by- 
level" analysis previously presented in Figures |22| |20| 



IX. RESULTS 

In Figure [24] we show the tttt i = elastic scattering 
phase shift for ttItt = 396 MeV as a function of center- 
of-momentum frame scattering momentum as extracted 
from finite-volume spectra. Discrete points correspond 
to a "level-by-level" analysis in which the Liischer equa- 
tion is solved for So{pcm) at each obtained Pcm with some 
justified assumptions made about the size of 62,4 at this 
scattering momentum, and with the degree of uncertainty 
about the higher £ partial waves reflected in a systematic 
error. The curves are the result of "global fits" to all the 
finite-volume energy levels assuming either an effective 
range parameterisation or just a scattering length, either 
of which are able to describe the energy spectrum well. 
The best estimates for the scattering length and effective 
range expressed in units of the pion mass on this lattice 
are 



but there is a very high degree of correlation (0.9) be- 
tween these values, and a pure scattering length of 
TOtt • ae^Q — —0.285 ± 0.006 can describe the data just 
as well. 

Figure [25] shows the tttt £ = 2 elastic scattering phase 
shift which can be well described by a scattering length of 
m^-ai^2 = (— 1.89±0.53)x 10~^. Statistically significant 
signals for elastic scattering in the £ — A wave were not 
observed and we estimate that m^ ■ |a^=4| < 1 x 10^^. 

We note here that the same L/us — 16, 20, 24 lattice 
ensembles (plus a larger L/ag — 32 ensemble) were used 
by NPLQCD to extract (5o(pcm) in [25 . They consid- 
ered many of the same frames, but limited themselves 
to the "scalar" irreps {Af for P = [0,0,0] and Ai for 
P ^ [0, 0, 0]), and they did not use a variational basis of 
operators. A comparison of results is shown in Figure [26] 
where low-lying levels are observed to have energies (and 
hence phase-shifts) that agree well, but where discrepan- 
cies appear at higher energies. The most significantly dis- 
crepant points (at {atPcm)"^ ~ 0.0017 and ~ 0.008) in the 
NPLQCD analysis correspond to levels which are either 
nearly degenerate with another level (the P = [0,1,1] 
ground statej^ or are highly excited {P = [0,0,0] sec- 
ond excited level). Since in our analysis we see no such 
discrepancies it may be that the variational method more 
reliably determines energies in cases where orthogonality 
of states is important. 

In our somewhat limited previous analysis of tttt I ~ 2 
scattering [^, we considered three pion masses and ob- 
served no significant dependence of the energy variation 
of 6q on the pion mass, which appeared to agree rather 



well with the experimental data. In Figure 27 we show 
our (5o,2(Pcm) obtained at m^ = 396 MeV along with the 
experimental data taken from [3]-(5]. Our data points 
have the absolute energy scale of the scattering momen- 
tum set using the J7-baryon mass procedure suggested in 

(SB, Pcm = {atPcu.){0^ with {atmn) = 0.2951(22) on 
these lattices [13]. Also shown is the tttt 7 = 2^ = 
phase-shift obtained using experimental information in 
multiple channels from a constrained analysis provided 
by the Roy equations, which implement manifestly cross- 
ing symmetry and the chiral behavior of the scattering 
amplitudes [34|l35]. 



X. SUMMARY 

A crucial step in the extraction of hadronic resonance 
properties is the determination of their resonant scatter- 
ing behavior. Within a Euclidean quantum field theory, 
the relevant elastic scattering matrix elements can be in- 
ferred indirectly through a systematic study of the spec- 
trum within a finite volume. In this paper, we extend our 



m^ -flfco = -0.307 ±0.013 
m^-rfco = -0.26 ±0.13, 



see Figure 14 



25 



0.002 



0.004 



0.006 0.008 



0.01 



0.012 



0.014 



o 



-10 



-20 



-30 



-40 



-50 




I 







P = [0, 0, 0] 


A L/a,, 


= 'H 






P = [0.0,1] 


° i/fls 


= 2U 


P= [0,1,1] 


°L/a, 


= 16 


p = [1.1, i; 






J 



(at ;?cm)^ 




scat. len. 
eff. range 



FIG. 24. Extracted 1 — 2 tttt elastic scattering phase-shift in S-wave, 5o{pcm), as obtained from analysis of finite- volume spectra 
with rriTy — 396 MeV. Center-of-momentum frame scattering momentum expressed in units of the temporal lattice spacing. The 
momentum region plotted is entirely elastic, with the 47r threshold opening at (atPcm)^ = 0.014. Colored points correspond to 
an analysis treating each energy level independently. The innermost errorbar is the statistical uncertainty, the middle errorbar 
refiects combined statistical uncertainty and uncertainty in {atm^, ^) and the outermost errorbar shows the total uncertainty 
including imperfect knowledge of ^2,4 (all errors added in quadrature). Curves indicate a global analysis of all energy levels 
describing the phase-shift by a scattering length or an effective range parameterisation. 



previous study [71 determining the i = and £ = 2 wave 
phase shifts in the tttt 1 — 2 systems, investigating more 
thoroughly the effects of operator basis, finite temporal 
extent, as well as the role of higher partial waves. 

With access to only modest lattice volumes, in order to 
map out the energy dependence with a significant num- 
ber of points, we determined the excited state spectrum 
in moving frames. This was achieved by constructing 
a basis of tttt operators transforming irreducibly under 
the reduced symmetry of a moving particle in a cubic 
box. Variational analysis of matrices of correlators built 
in this operator basis leads to extraction of excited state 
energies with high precision. 

The construction of a basis of operators with suit- 
able quantum numbers corresponding to the tttt system 
in-flight is a significant extension beyond the previous 
work, and has allowed for the determination of the phase 
shifts at many discrete energies below the 47r inelastic 
threshold. This increased operator basis, covering many 
irreducible representations, allows for more constraints 
on the contributions of higher partial waves. However, 
the weakness of tttt scattering in the isospin-2 channel 
presents a particular challenge to extraction from finite- 
volume methods. The changes in energy with respect to 
non-interacting pions determine the phase-shift and since 
these are small it is important to take care over system- 
atic effects that may be small in absolute terms but which 



could be large on the scale of the energy shifts. 

We reduced the contribution of excited pion-like me- 
son states to our tttt correlators by using optimised pion 
operators. These operators are constructed from a lin- 
ear combination of composite QCD operators with pion 
quantum numbers and their important property is that 
they relax to the ground state faster than any single sim- 
ple operator construction. The reduced contribution of 
tttt"^ states to our correlators allows analysis at earlier 
Euclidean times. 

At larger Euclidean times, the effect of the finite tem- 
poral extent of the lattice can be observed, distorting the 
time-dependence from the desired sum of exponentials 
corresponding to discrete state energies. We have explic- 
itly accounted for the largest unwanted finite-T effects 
leaving sub-leading effects which are somewhat smaller 
than the statistical uncertainty. 

The reduced symmetry of a cubic box at rest is such 
that 60 always appears with some sensitivity to ^4, but 
the very small value of 64 throughout the elastic region is 
such that the rest-frame spectrum is mostly independent 
of ^4. On the other hand, the symmetry of a cubic box is 
further reduced when placed in-flight and So extractions 
become sensitive to the value of 62, which is not neces- 
sarily negligibly small. We investigated the effects that 
non-zero values of S2.4 can have on the finite- volume spec- 
trum using a toy model showing that some energy levels 
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FIG. 25. Extracted 7 = 2 tttt elastic scattering phase-shift in D-wave, 52(pcm), as obtained from analysis of finite-volume 
spectra with m,r = 396]VIeV. Center-of-momentum frame scattering momentum expressed in units of the temporal lattice 
spacing. Momentum region plotted is entirely elastic, with the A-k threshold opening at (atPcm)^ ~ 0.014. Colored points 
correspond to an analysis treating energy regions locally as described earlier in the manuscript. The inner errorbar is the 
statistical uncertainty, and the outer errorbar reflects the combined statistical uncertainty and uncertainty in atm^, ^ and the 
value of ^4 (errors added in quadrature). Curves indicate a global analysis of all energy levels describing the phase-shift by a 
scattering length. 
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FIG. 26. Our So extraction (colored points) compared with those of NPLQCD (grey points) over the elastic region (left) and 
zoomed in to small scattering momentum (right). 



can show significant sensitivity. 

We attempted to account for the effects of higher par- 
tial waves on the extraction of <5o,2j finding that they are 
generally small (except in a limited number of sensitive 
cases identified in the toy model analysis) . We associated 
a systematic error with our imperfect knowledge of them 
that was found to be always smaller than the statistical 
uncertainty. We found that the finite volume energies 



could be well described by a scattering length parame- 
terization in both i = and £ = 2 over the elastic region. 
The fit could be moderately improved by adding an effec- 
tive range in £ = 0, albeit with a significant correlation 
between the effective range and scattering length. The 
fits did not indicate the need for significant strength in 
the £ = 4 wave. 

The calculations reported in this paper were performed 
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FIG. 27. Extracted 1 = 2 tttt elastic S'-wave(red), Z)-wave(blue) scattering phase shift (for m^ = 396 MeV, all errors combined). 
Shown in grey the experimental data from [SHS] and the constrained analysis using Roy equations [341 1351 (black line, grey 
band). For the heavy pion mass the entire region plotted is elastic while for the experimental data only p^^ < 0.058 GeV is 
elastic. 



at only a single pion mass of 396 MeV. While they demon- 
strate that the procedure outlined can indeed determine 
scattering phase shifts with a high degree of confidence, 
the obtained results cannot be directly compared with ex- 
perimental data. Future calculations using lighter pion 
masses will be required, as will eventual consideration of 
other systematic effects such as the lattice spacing de- 
pendence. The results presented in this paper supersede 
those presented in [7] which considered only rest-frame 
correlators using un-optimised pion operators and where 
finite-T effects were not fully accounted for. 

The techniques developed in this calculation are a nec- 
essary ingredient to future investigations of resonances 
in hadron-hadron scattering that arise from the strong 
interactions. At unphysical pion masses, the phase space 
available for decays can be small as seen in studies of the 
/ = 1 TTTT sector [36H39] giving rise to a rapid variation 
of phase-shift with energy. Thus, the formalism and con- 
struction of operators in-flight developed in this work will 
be necessary to compute a sufficient number of energies 
within the resonance region to allow for a reliable deter- 
mination of resonance parameters. To compute these en- 
ergies, the operator basis used in the variational method 
will feature both single and multi-hadron constructions. 



Annihilation diagrams will arise, which as shown in the 
isoscalar meson sector p^ , can be efficiently constructed 
using the "distillation" method. 
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Appendix A: Multi-particle operators 

In this appendix we give a generalisation of and fur- 
ther details for our multi-particle operator construction. 
In Section [A l| we describe our general multi-particle op- 
erator construction and in Section lX^ discuss the method 
used to calculate the Clebsch-Gordan coefficients. In 
Section IA3I we state our conventions for the lattice ir- 
reps and choices of rotations before giving some exam- 
ple Clebsch-Gordan coefficients in Section |A4[ Further 
Clebsch-Gordan coefficients are given in supplementary 
material 1431. 



tion coefficients, Sj\, which were discussed in Ref. [27] . 
are given in Table IX Note that here we give subductions 
coefficients for creation helicity operators which trans- 
form like states, i.e. under a rotation by (j) about the 
axis defined by p they transform as R^ |A) — e^"^'^ |A). 
In Ref. [22] we showed that these subduced helicity op- 
erators are useful for studying mesons with non-zero mo- 
mentum on the lattice. 

In general, a two-particle creation operator with total 
momentum P can be constructed from the product of 
two single-particle operators, 



1. Operator construction 

Here we describe our construction of multi-particle op- 
erators, generalising the discussion in Section [ll] We will 
consider two-particle operators but the procedure can be 
applied iteratively to construct multi-particle operators 
for a larger number of particles. 

In Refs. (mill] we discussed how operators with a def- 
inite continuum J^ and spin-component M, O'^ •^'\p = 
0), can be constructed out of gauge- covariant derivatives 
and Dirac gamma matrices. The appropriate lattice op- 
erators were formed by subducing these continuum oper- 
ators into octahedral group irreps, 



0^1(0) 



-E^h 



o-^ (0) 



t 



M 



where A^ is an irrep of 0/j , /i is the irrep row and S j j(^ 
are subduction coefficients discussed in Refs. [TTlfT^ITl] . 
In Ref. j27j we discussed how to construct helicity op- 
erators^ 



0-' ^^{p) 
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'Di^m 



0-' 



\v) 
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(Al) 



where A is the helicity, J^ refers to the spin and parity 
of the operator with p = and Z? is a Wigner-D matrix; 
a summary of our conventions is given in Appendix |A3| 
and we refer to Ref. [TT for more details, for example, 
on the choice of i?. From these we constructed subduced 
helicity operators, 
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A.fi 



O' '^(p) 



where A is an irrep of LG(p) (the little group for momen- 
tum p), fi is the irrep row and fj = P(— 1)"'. The subdue- 
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TABLE IX. Subduction coefficients, 5^' , for integer spin. 



f,,A 



jAj < 4 and with s ■ 
text. 



sign(A); other notation is defined in the 



here we give expressions for creation operators 
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j[Ai{fcir;A2{fe2r] 



(P) 



E 

fcie{fci}* 

k2e{k2y 

ki+k2=P 



C(PA/^;fciAi/ii;fc2A2/i2) Oai,.i(A:i) OA2^.Ak2) 



where Ai_2(mi,2) and A(/l(,) are respectively irreps(irrep rows) of LG(A:i^2) and LG(P). C are the Clebsch-Gordan 
coefHcients for Ai({fci}*) (g) A2({fc2}*) — > ^{P) which we discuss in the following section. The sum over ki2 is a sum 
over all momenta in the stars of A:i_2, {^1,2}*, i-e. all momenta related to fci.2 by an allowed lattice rotation. In other 
words, the sum is over i?fci.2 V i? G O^ ; the restriction that fci + fc2 = ^ is equivalent to requiring R e LG(P). 



r 



2. Induced representation method for calculating 
Clebsch-Gordan coefficients 

We use the projection formula with the induced rep- 
resentation [35 to construct the Clebsch-Gordan coeffi- 
cients for Ai (g) A2 ^- A where Ai, A2 and A are each 
irreps of, respectively, groups Gi, G2 and G. For our 
purposes, these groups will be the double cover of the 
octahedral group, O^, or a httle group C O^^. 

A group, G, can be partitioned into cosets by a sub- 
group H C G. Two elements x, y € G are in the same left 
coseipj if and only if y~^x G H. One coset, containing 
the identity element, will be the subgroup H itself. The 
number of cosets n = |G|/|H|. A coset representative is 
one element from the coset; Ri,R2, . . . ,Rn are a set of n 
coset representatives, one from each coset (a transversal). 

If A is a I A|-dimensional irrep of H, then T = A(II) f G 
is a (n|A|)-dimensional unitary representation of G in- 
duced from irrep A of H. It is defined for i? € G by 



-L [-tt/irjs — 







RRj)j- 



if i?j RRj 
otherwise 



eH 



(A2) 



Here i,j label the coset and r, s the rows and columns of 
the irrep A. 

To make this more concrete, consider the double-cover 
octahedral group with parity, G — O^ (|G| = 96), and 



the little group for momentum p^cf = [0,0, 1], H = Dic4 
(|II| = 16). There must be 6 left cosets with this lit- 
tle group and these correspond to the 6 momenta in 
{Prci}*- If ^1 and i?2 are in the same coset, then 

R2Pxci — RlRi R^Pvci = RlRnPrei whcrC Rh & H. But 

from the definition of the little group, RhPtcI — Prcf, 
and therefore i?2Prcf — RiPicf- The converse can also 
be shown to be true. Therefore, each left coset can 
be labelled by p G {prcf}* with elements R such that 
Rpid = P- 

Therefore, in Eq. |A2[ the indices i,j refer to a par- 
ticular momentum direction, p g {pref}* ■ In effect the 
induced representation splits up R into a piece Rh in 
the little group [giving A{Rh)] and a piece which rotates 
the momentum direction from p^^i to p. The freedom 
to choose a particular coset representative for each coset 
is the same freedom as the choice of lattice rotation R 
which rotates Prct to p (see Ref. [21] )• It is not important 
which particular element is chosen as the representative 
but this choice should be made consistently; we discuss 
our conventions in Section [A 31 

Once the induced representations, F^ — Aj(Gi) t O^ 
{i — 1,2), have been constructed, the Clebsch-Gordan 
coefficients can be generated using the projection formula 
in the same way as Clebsch-Gordan coefficients for O^ 
0/j ^ 0/j . The projection formula gives 
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[Oa,,;(M*J]^[Oa,,,(M,\)]^ 



(A3) 



where P = {pi}^^ + {P2}j2 is fixed, fi, p,^P , p.2 label ir- 
rep rows, and Ji,ji and J2,j2 label the different mo- 
menta in the sets {pi}* and {^2}* respectively. Note 



r 



that this is written for operators which transform as 
RO\^ = T,t,'HR)t^'t^'^A^'- After forming the appro- 
priate linearly independent combinationfr^ the Clebsch- 



^"^ note that if H is not a normal subgroup, the left and right cosets 
are different 



^^ linear combinations of ^',^1,^2', it is sufficient to consider one 
particular ii and 12 
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Gordan coefficients can be read off up to a phase choice 
and normalisation. We choose the phase so that our re- 
sulting correlators are real. 



When 



P - Wl + Wl - 0, G 






and the 



sum is over all R 
consider 11=^2 = 
coefficients. 



e 0,j . In this case it is sufficient to 
1 to generate all the Glebsch- Gordan 



When P — {pi}*^ + {P2}|2 7^ 0, G is a little group and 
the sum over R is restricted to those elements in the little 
group. A particular choice of «i and 12 can be made such 
that P — {pi}*^ + {p2}*2 ; because R is in the little group, 
the rotations of pi and p2 will automatically ensure that 
P is fixed. 



Lattice and little group conventions and 
rotations 



As described in Rcf. '27J , we break the (active) rotation 
R appearing in Eq. |A1| into two stages: R = i?iat-Rrcf- 
The first rotation, R^eU takes (0,0, \p\) into Prof, where 
Prcf is a reference direction for momenta of type p (i.e. for 
{p}*)- In Table Ix] we give the specific rotations that we 
use for i?,.cf. We use the same convention as in Ref. [27], 
namely that a rotation -R^,^,^ = e~^'^'-'^ e~^^'^y e~^^'-''' ro- 
tates around the z-axis by ipi then around the y-axis by 
9 and finally around the z-axis by (j> (with a fixed coordi- 
nate system). In Table XI we give the rotations, i?iat, for 



each momentum of the form [0, 0, n]; these correspond to 
the coset representatives discussed in Section |A2| Ro- 
tations, i?iat, for other types of momenta are given in 
supplementary material |43j . 

In Tables pCTTl pClTIl pClVl and pTV] we give our choice of 
representation matrices for, respectively, the little groups 
Dic4, Dic2, Dic3 and C4. For Dic2, Dica and G4, the 
rotations and reflections refer to a coordinate system 
which has been transformed using -Ricf, so that p de- 
fines the new z-axis. Note that the convention is such 
that states in little group irrep A (row /i) transform as 

^|AA^>=EM'A^,JAAi'). 
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TABLE X. Rotations, R^ci, used, as described in the text; 
liere n is a non-zero integer. 

4. Clebsch-Gordan coefficients 

Clebsch-Gordan coefficients for zero total momentum 
[P — ki + k2 = 0) with fci = A:2 = are given in 
Ref. [28] . Example Clebsch-Gordan coefficie nts for P = 

[xvi|xviii 



with ki = — ^2 7^ are presented in Tables 
others are given in supplementary material 143J . We show 
Clebsch-Gordan coefficients for Ai ® Ai ^ A; those for 
A2 ® A2 are identical, and those for Ai ® A2 and A2 ® Ai 
follow by switching the target irrep's parity, A='= — >■ A^. 

The cases of non-zero total momentum where fci = 
and P = k2 ^ {oT k2 = and P = fci ^ 0) are trivial 
because there is only one momentum in the sum. We 
have 



C ( A(^) [P] ® A+ [0, 0,0]^ A{^l) [P] 



1 



where /z is the irrep row; those for other target irreps and 
rows are zero. In addition, 

C f A(Ai) [P] A^ [0, 0, 0] ^ A' i^l')[P]) = 1 , 



where if A = Ai,A2,Bi,B2 then A' = A2,Ai,B2,Bi 
respectively and if A(/^) = E2{l),E2i2) then A'(/i') = 
-£^2(2), £^2(1) respectively (the rows are swapped around). 
The coefficients for other target irreps and rows are zero. 
For non-zero total momentum and ki,k2 ^ we 
present some examples of Clebsch-Gordan coefficients 
in Table [XIX[ others are given in supplementary ma- 
terial 33^. We show Clebsch-Gordan coefficients for 
A\ ® A\ — >■ A; those for A2 ® A2 are identical, and 
those for A-y ® A2 and A2 ® A-y follow by replacing 
the target irrep A(^) = Ai, y42,Pi, P2,^2(l),£'2(2) by 
A(^) = A2,^i,P2,Si,^2(2),£;2(l) respectively. 
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TABLE XL Rotations, -Riat, used for momenta of type [0,0,n], as described in the text; here n is a non-zero integer. 
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TABLE XIL Choice of representation matrices for the Dic4 Httle group. / denotes the identify transformation, R{(t)) denotes a 
rotation around the z-eads by and 11 denotes a reflection in the yz plane (x — >■ —x). Note that, because we are considering 
only irreps relevant for integer spin, the representation matrices for R{(f) + 2ii) are the same as those for R((f>). 
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TABLE XIIL As Table [Xn] but for the Dic2 little group. 
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TABLE XIV. As Table IXlT] but for the Dies little group. 
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TABLE XV. As Table IXlT] but for the C4 little group. 
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TABLE XVI. Clebsch-Gordan coefhcients for yli[0,0,n] ® 
j4i[0,0, n] — >■ A[0,0,0] where /i is the row of A and n is a 
non-zero integer. 
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TABLE XVII. Clebsch-Gordan coefficients for Ai[0,n,n] ® Ai[0,n,n] 
integer. 



A[0, 0, 0] where /i is the row of A and n is a non-zero 
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TABLE XVIII. Clebsch-Gordan coefficients for A\ [n, n, n] ® Ai [n, n, n] 
integer. 



A[0, 0, 0] where /i is the row of A and n is a non-zero 
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TABLE XIX. Clebsch-Gordan coefficients for Ai[0, -n,0] ® ^i[0,n,n] 
A[0, 0, n] (right) where p is the row of A and n is a non-zero integer. 



A[0,0,n] (left) and Ai[-n, -n,0] ® Ai[n, n, n] 



